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<N ; I. INTRODUCTION 
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^ ■ In this paper we use the phenomenological model of Ref. [jj (KKT model) to analyze 

j> ! inclusive hadron production in pA collisions at RHIC and LHC. We focus on energy, rapidity, 
^ ' centrality and transverse momentum dependence of the nuclear modification factor Rp^ {p 
. . . stands for either proton or deuteron) for vr, D and B meson production. Theoretical analysis 

reveals that coherent effects in hadron production become phenomenologically 
significant at RHIC. They consists of (in the nucleus rest frame): (i) multiple rescattering of 
the incoming parton system in the nucleus. This is controlled by the parameter a'^A^^^ ~ 1. 
And (ii) small-x quantum evolution controlled by the parameter = as ln(l/x) ~ 1. Both 
effects lead to the gluon saturation of a dense parton system at small x which manifests itself 
in a different ways depending on the kinematic region we are interested in: (i) at not too 
small X (e~^/°^ < x <^ 1) it contributes to the Cronin effect as the result of multiple 
rescatterings of the projectile in the target; (ii) at a; < e"^/"'' the small- a; evolution starts 
off which results in a dramatic suppression of at forward rapidities at RHIC. Similar 



effects have been predicted 



10[ | and observed ll| for open charm production. It is 



of great interest to observe the evolution of these phenomena with energy. Thus, motivated 
by the approaching start of the nuclear program at LHC we decided to combine calculations 
of inclusive hadron production in a framework of a single model. 

The paper is structured as follows. In Sec. [Ill we review a phenomenological model for 
the dipole scattering amplitude N{r, y) [l| which has been previously used to successfully 
describe the data on light hadron production at RHIC. Althou gh N(r, y) can in principle 



be analytically calculated using the small-x evolution equation [1^, Il3| , its a very difficult 
and not yet solved problem. In addition, there is a fare amount of uncertainty in how the 
NLO terms in small-x evolution change the behavior of N{r, y) as well as in the effect of 
fluctuations in the dense partonic system. On the other hand, rather general arguments 
allow identification of a few key features which hold for N{r,y) in a dense parton system. 
This features are heeded in constructing a phenomenological models for N{r,y). In this 
perspective we argue in favor of the KKT model of Ref. Sj. 

In Sec. mil and Sec. IIVI we present theoretical results which together with the model of 
Sec. [Tl]enables calculation of pion and open charm and beauty production at RHIC and LHC. 
These results are summarized in figures Fig. [H Fig. [3] and Fig. |H They teach us that at LHC 
at rapidities y > 0, the nuclear modification factor is a very slow function of rapidity. This is 
an anticipated result. Indeed, the amount of suppression can be estimated as RpA — 1/N^~i^, 
where 7 is a depends on rapidity and transverse momentum 2l, ISl]. Eqs. ([2]) and (flTl) imply 
that 7 decreases with energy and rapidity and increases with transverse momentum and 
mass, although it is a very slow function of its variables. Thus, we observe almost the same 
suppression pattern as the function of transverse momentum, rapidity and centrality for 
gluons and heavy quarks. To emphasize the mass and transverse momentum dependence 
of RpA we plotted it versus m in Fig. [5] up to the top-quark mass. Despite that x is small 
in the kinematic region of this figure, the geometric scaling and hence the KKT model, are 
expected to break down. Still there is a vague theoretical understanding of how exactly it 
occurs. This too is discussed in Sec. HVl 
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II. A MODEL 



In the dipole model jl4|, cross sections for inclusive hadron production can be expressed 
through the gluon dipole forward scattering amplitude N(r,y), where r is the transverse 
separation and y is rapidity. This quantity can be calculated for any r andj/ > 1 using the 
small-x evolution equations, the most useful of which is the BK equation fl2l. \l3 \ . In order 



15|, 



16 



18 



19 



one needs to 



to be able to utilize an exact numerical solution to BK 
(i) include the NLO corrections into the BK equation and (ii) perform global analysis of all 
small X data. Until this program is carried out, the most practical way to proceed is to use 

solution to the BK equation. 



a model for N(r,y) which has the most essential properties o: 



Several such models have been suggested l21|, l22|, l23|, |2j, |25|, |26|. Although they are 
not very much different from one another, the recent analysis of the numerical solution to 
BK equation performed in Ref. 27| seems to imply that the KKT-like models are closer to 
the numerical solution than the other models. The main observation of Ref. [27] is that the 
anomalous dimension 7 at the saturation scale is about 0.44. This should be compared to the 
anomalous dimension of 0.5 used by models based on the double logarithmic approximation, 
such as KKT, and to the anomalous dimension of 0.628 used by models based on the saddle- 
point approximation. It should be kept in mind that 7 is the most important parameter 
determining the maximal possible suppression of the nuclear modification factor due to the 
gluon saturation effect. On the other hand, the KKT model is certainly oversimplified and, 
as more accurate data appear, it will be replaced by a more realistic models. 

For the reasons described above and encouraged by success of the KKT model in describ- 
ing inclusive light hadron production at RHIC we set to employ it for calculation of inclusive 
hadron production at LHC According to Ref. {l| the forward gluon dipole scattering am- 
plitude is parameterized as follows 



^{.11: y) = ^ ~ exp 



(1) 



where we denoted \r\ = r. The anomalous dimension is parameterized in such a way as to 
satisfy the analytically well-known limits of (i) r — > 0, y fixed and (ii) y — > 00, r fixed: 



l{r,y) 



1(1 + 



|5(r,,j)| + V2|5(r,3/)|+28C(3)^' ^ - ' 

1 y <yo, 



(2) 
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where 

In the double logarithmic approximation we can replace ^.Xj (4A;^). The gluon saturation 
scale is given by 

g2(y) = A2 ^1/3 gAj; _ g^^3 Qgy2 ^Aj; ^^^^^ ^ 



Parameters A = 0.6 GeV and A = 0.3 are fixed by DIS data 2l|]. The minimal saturation 



scale used in ([3]) is defined by Q'j.Q = Qliyo) with the value of rapidity at which the small-x 



quantum evolution effects set in. Fit to the RHIC data yields uq = 0.5 [l] ^. 

In the quasi-classical approximation (at 7 = 1) the forward scattering amplitude of the 
quark dipole A^Q(r, is specified by the same formula with the gluon saturation scale 
rescaled by the color factor Cp/Nc- Therefore, we model NQ{r_,y) in the whole kinematic 
region by 

NQ{r,y) = 1 -exp |-l(r2g2/2)7(n.)| , (5) 
where in the large Nc limit Cp/Nc = 1/2. 

III. LIGHT HADRON PRODUCTION 

The light hadron production cross section is a sum of two terms: inclusive gluon and 
valence quark production and hadronization. Inclusive gluon production dominates in most 
of the kinematic region safe for the most forward rapidities in the proton/deuteron frag- 
mentation region, where Bjorken x of nucleus [xa) acquires its lowest possible value for a 
given ^/s, while that of proton/deuteron {xp) is close to unity. In that region rescattering 
of valence quarks of proton in nucleus dominate to the hadron production cross section. 

The cross section for inclusive gluon production is given by 

_^g^oi^fSa x(^^_ Y I dzTJoikTZr) \n^d^^[zTd;,^NG{zT,y)], (6) 
d^kdy TT^ ^ Jq 2^/1 

where is a scale associated with deuteron and is fixed at fi = 1 GeV thereof. Expression 
Xp^{l — Xp)'^ is a model for the gluon pdf of proton. Inclusive valence quark production cross 

^ In [l| it has been argued that a non-perturbative scale ( "intrinsic fcr" ) which is responsible for the Cronin 
effect at low energies has a little importance at RHIC. It is therefore neglected in this paper. 
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section [28 1 

The nuclear modification factor is defined as 

da^(pA) dN{pA) 

R^^(k, y) = = ''""fl (8) 

^ d^kdy ^^c'^^^d^kdy 

where '^^^^^^y and '''^jf^y are multiphcities of hadrons per unit of phase space in pA and pp 
collisions. Both expressions for gluon ([6]) and valence quark ([7j) production contribute to 
the hadron production cross section in pA and pp collisions. The cross section of hadron 
production reads 

da^ f dz dac ^ , , 



d^kdy J z'^ d^kdy 
dz dijQ 
z'^ (P'k 



{k/z) xqviy, k/z) Dq{z, kr) F{k/z, y). (9) 



We use the LO fragmentation functions from Ref. j29(] with the renormalization scale of the 
fragmentation functions equal k. 

Equation ([7j) is derived for production of a valence quark in the proton/deuteron frag- 
mentation region. To generalize it to smaller values of Bjorken x one has to convolute it 
with the proton's/deuteron's valence quark distribution, which is fixed by quark counting 
rules at high x and by the leading Regge trajectory at low x 

xqv{x) = l.m{l-Xpfxl\ (10) 

I — I 

where Xp = {kx/ \/s) e^. Valence quarks are increasingly less important at low x [36], where 
the quark production is dominated by gluons splitting in qq pairs. The factor of x^^ insures 
that this is indeed the case here 36|. Analogously, the high x behavior of the nuclear gluon 



distribution is taken into account by introducing the function F{k,y) 

Fik,y) = il-x^r —--] . (11) 



k^ + A\ 

where the Bjorken x of a gluon in the nuclear wave function is given by xa = {k/ ^) 
and as = 0.3. The last factor in Eq. ( fTTi) arises when we impose momentum conservation 
constraint on the anomalous dimension of the distribution functions. 

Eq. (fTTj) is an empiric way to include the higher order pQCD corrections. It is well-known 
that the NLO pQCD corrections are required to describe the inclusive 7r° production in pp 



6 



collisions at RHIC 



30 



3l|. Therefore, these corrections must be taken into account both 



in the numerator and the denominator of Eq. ([8]). Unfortunately, the theoretical work on 
inclusion of the NLO corrections in the gluon saturation regime proved to be a very laborious 



problem. Only recently, calculations have been done 



corrections into the forward scattering amplitude 32 



or inclusion of some running coupling 



33 



351]. This is the main reason 



why presently, we have to resort to the phenomenological approaches to the NLO corrections 
such as the one described here. 

The results of numerical calculations are exhibited in Fig. [H As has been already pointed 




FIG. 1: Nuclear modification factor for pion production at RHIC and LHC. 

out in pj], since rapidity r] = 3.2 at RHIC corresponds to almost the same x as ?7 = at LHC, 
the nuclear modification factor is expected to be the same. This seen in Fig. [H We also 
observe that the difference between rj = and rj = 6 is less than 20%. This teaches us that 
fast onset of the gluon saturation observed at RHIC is replaced by virtually y-independent 
behavior at LHC (from r] = to r] = 6). Behavior of RpA at rapidity r/LHC = —3.3 and 
rapidity ?7rhic = is expected to be similar. The amount of suppression perceived in Fig. [1] 
depends mostly on centrality, i.e. Ai'cou: RpA — A~^l^ ^ ^coiV^- This yields R^a = 0.60 for 
minimal bias events and RpA = 0.46 for 0-20% centrality cut in agreement with Fig. [H 

To conclude this section we would like to emphasize that we explicitly neglect a possible 
effect of gluon saturation in a proton or deuteron. This approximation is perhaps sufficiently 
good for the nuclear modification factor which involves the ratio of the cross sections and is 
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intentionally constructed to look for nuclear effects. 



IV. HEAVY QUARK PRODUCTION 

Production of heavy quarks at small x is also affected by gluon saturation in a way similar 
to that of gluons [l^. The main difference, however, is that the effect of gluon saturation 
is postponed to higher energies/rapidities for heavier quarks as compared to lighter quarks 
and gluons. This is because the relevant x is proportional to m_|_ ~ (m^ + k'j_Y^'^ and hence 
is higher for heavier quarks at the same values of ^/s, y,k. 

In this section, we first of all review the theoretical approach to the qq pair production 



at small x, see Ref. 



9|. Let us introduce the following notations: k is the produced quark 
transverse momentum, q is the gluon transverse momentum, a = is a fraction of the 

light-cone momentum of gluon carried by the produced quark; x^ and are the transverse 
coordinates of the produced quark in the amplitude and in the complex conjugated amplitude 
respectively; and y^ are the corresponding coordinates of the antiquark. Transverse 
coordinates of gluon in the amplitude u and in the complex conjugated amplitude v are 
given by M = a Xi + {1 — a) X2 with u = \u\ and x 12 = Xi—X2, {xu = \xi2\) and analogously 
'or V. With these notations the single inclusive quark production cross section is given by 







da 



^—7 / (fxi(fx2(fyid'^y2S{u-v-a{xi2-y,n)) [ c^^e '-^-^ (12) 
'^T^r J Jo 



d'^kdyd% 2(2 7r 

3 

X ^ ^ ^ ij i^^l-i ■L2-I y_i^ ^2-1 '^ij^-!^^ (13) 

where the delta function comes about after integration over the antiquark transverse mo- 
mentum which implies x_2 =1/2' '^^^ products of the light-cone "wave functions" are detailed 
as follows [91] see Fig. [2] 

^ii{.x^.X2;y ,y -.a) =ACf(—] I ^^2(^1,^2! «) ^2(y,,l/„; a) ^ [{l-2af 

X (^12 ■^) (^12 + (eii«i^i2i) {^kiVkyi2i)] + Fi{x^,X2;a) Fi{y^,y^;a)m 



u ■ V 



+4:a'^ (1 - af Fo{x^,X2] a) Fo{y^,y^; a) - 2 a (1 - a) (1 - 2 a) 



X12 u 
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3> 
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22 





33 



3> 





12 





23 



3> 




13 



3> 




FIG. 2: The diagrams contributing to the quark-anti-quark pair production in the quasi-classical 
approximation. Disconnected t-channel gluon lines imply summation over all possible connections 
to the adjacent s-channel quark and gluon lines. Adopted from Q]. 



V ■ V 
^-12 - 



yuv 



(14a) 



$22 (xi,X2;?/^, =4Cf( — ) \ Kiimxu) Ki^myu) [(1 -2af 

-1-2 V ^ / I ^12 2/12 



11 • 1) I 

X {xi2 -u) {y.o + {(^ij 'U'i X 12 j) (ekiVkVui)] + Ko{mxi2) Ko{myi2) > , (14b) 

^i2ix^,X2;y^,y-,a) = -ACp (—] "^ |i^2fe,^2; «) ^i("^l/i2) ^ ^ [(1 -^af 

-1 \TX J y Xi2 2/12 U 

11 . 1) 

X fe2 (2/1, - I!) + {(^ijUiXuj) (ekiVkVui)] + m Fi{x^, X2; a) Koimyu) 

—IZ y yZ 

y ■ y_ ^ 

-2a(l - a) (1 -2a)=i^Fo(xi,X2;a)Ki(m2/i2) ^ , (14c) 

$33(^1, ^2; ^1)^2'") " "^11(^1,^2; Mi>^2'") + '^22(^1, ^2; ^1)^2'") + '^12(^1, ^2; ^2'") 

+ $2i(xi,X2;y^,y2;") (14d) 
$13(^1, ^2; ^2'") = -"^11(^1, ^2;^i>y2' ~ "^12(^1, X2;y^,y2; a) (14e) 

*23(a,^2;yi>^2'") = -'^2l(^l,^2;^i,^2'") ~ '^22(^l,^2;^i,^2'") *^^^^) 

^ij{xi,X2]y_^,y_^,Q) = ^ji{y^,y^;x^,X2;a). (14g) 
The auxiliary functions Fi, F2 and Fq are defined as 

^2(0;]^, X2; a) = J dq Ji{qu) Ki (^12 a/ w? + a{l — a)^ a/ + a{l — a) , (15a) 
-^1(^15^2;") = y dqJi{qu) Kq(^i2 v^m^ + g2Q,(^i _ ^ (^15b) 
Fo{x^,X2;a) = dqq jQ^qu) Kq(^12 ^/ m? + q^ a{l - a)^ , (15c) 

where m = a;]^2 = iZIi ~ 2:2, a;i2 = |aii2l! and q = + ^2- 

Using the definition of the gluon saturation scale Qs (see (|4])) 

Ql = 47ralpT{b) (16) 

with p the nucleon number density in the nucleus and T{b) the nuclear profile function, we 
write (in the large Nc approximation) 

^iiix^,X2;y^,y^;a) = e-^(^i-^i)''3^n(i/|x,-?;j^)-i{^2-i/2)'Q^n(i/ix2-2^,l/^) ^ (17a) 

^22{x„X2;y^,y^;a) = e-i(^-^)'Q^ ^'^(i/l^"^!^) , (17b) 

'^33ix^,X2;y^,y^;a) = 1, (17c) 

^i2{x„X2;y^,y^;a) = e-i(^i-^)'«' Mi/|x,-.1m) ^ (lY^i) 

^23ix„X2;y^,y^;a) = e-i"^Q^-(V-A.) , (17e) 

^i3ix^,X2;y^,y^;a) = e-i-i«' in(i/a=2^) (^i^i) 
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All other Sj^ 's can be found from the components listed in (fT7|) using 

'^ij{xi,X2;y^,y^;a) = 'B.ji{y^,y^;x^,X2;a) (18) 

similar to (|14gp. 



If the typical gluon momentum q is much smaller than the produced quark mass, the 
above expressions can be significantly simplified. Indeed, since a{l — a) < 1/4 we get 

(fa{l - a) < , (19) 

and the auxiliary functions read 

F2{xi,X2;a) = Ki{xi2m) mu^^ , (20a) 
Fi{xi,X2'^a) = KQ{xi2m) u'-^ , (20b) 
Fo{x„X2;a) = 0. (20c) 

In this approximation the only non-vanishing products of "wave functions" are given by 

^nixi,X2;y^,y^;a) = ^22ix^,X2;y^,y^;a) = -^i2ix^,X2;y^,y^;a) 



ACf ( — ) \Ki{xi2m)Ki{yi2m) [(1 -2af {x^2 ' m) (z/io ' v) 



XL ' V 

+ {eijUiXi2j) {tki Vk yi2i)] + Ko{xi2m) Koiyum^ 



y2 y2 



Averaging over all directions of gluon emission from the valence quark using 
{eij Ui Xi2j ekiVkVui) = (1/2) uvxi2 ■ y^^ we arrive at the well-known result j3, Q 

^ii{x^,X2]y_^,y_^,a) = 













1 uv 1 



ACf{ — \ —^±^^^[{l-a') + a']K,{xum)Ki{yi2m)+Ko{xum)Ko{yi2m)\ (21) 

Approximation (1191) . employed to derive fl^Tl) . amounts to the /cT-factorization of the 
gluon distribution of proton (see fl2^ ). It breaks down if x is so small that proton can no 
longer be regarded as a dilute object: it acquires a hard scale ~ which increases with 
energy. However, in such a case even a more general formulas (IT^ . (ITTI) which go beyond 
the /cT-factorization of proton by do not include the high density effects, become invalid too 
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as they neglect the nonlinear gluon evolution in proton. Discussion of the gluon saturation 
in proton is beyond the scope of this paper. 

It is convenient to express the scattering amphtude in terms of the vectors Xi2, y^^, u 
and y_ and integrate in f[T^ over these variables. We have 

^i-y_^ = ^12 - y_i2 ' (22a) 

a-^i = 0, (22b) 

Xi-v = Xi2-a^^2' (22c) 

= -«yi2' (22d) 

y^-u = -a + ' (22e) 

y^-u = -ax^2- (22f) 

Thus, sum over all rescattering factors including the signs of $ij's is given by (we omit 
logarithms ln(l/x/i) for brevity) 

E{x,2,y,2,a) = e-i(^i2-^,,)^{Q?/2) +e-|«^fe.-I,,)^(Q?/2) 

Using the delta function in (fT3|) to integrate over u and integrating over y_ in the leading 
logarithmic approximation we derive the final result 

da _ CFaW [ ,2^^,2„^^^-^k■(.,,-y^ 



(Pkdy(Ph 47r5 



j^daj rf2xi2Ai2e-*^ (^--^i2) ln(l//x|x,2-y^J) 



=i^-=i^[(l - a^) + a^] Ki{xi2m) Ki{yi2m) + Ko{xi2m)Ko{yi2m) \ E{x^2^ y ; a) (24) 

Before we turn to the numerical results, we would like to remark on the large Nc cor- 
rections to ( l24l) . In our approach these corrections arise in the scattering amplitudes ( fTTIl . 
!5xact expressions for those amplitudes in the quasi-classical approximation were derived in 
3, S]. Numerical uncertainty due to omission of the ~ terms is expected to be of 



the order of 10%. A detailed analysis performed in 



suggests that these corrections are 



even smaller. Another source of the large corrections is the high energy gluon evolution 



beyond the mean-field approximation. These corrections were studied in detail in [40l] and 
are shown to be at the 1% level. Alluding to the poorly known NLO corrections as well 
as the accuracy of the present experimental data, we think that neglecting the large 
corrections is justified. 
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Numerical calculations are performed along the same steps as in Sec. IIIII we replace the 
dipole scattering amplitude by the model ([1]) and ([5]) and take into account the kinematic 
factors ( ITOl) . (|TT1) . In ([2]) and ([3]) we substitute ^ l/{4m\). To obtain spectra of D and B 
mesons we convolute the cross section ( 12^ with the Peterson fragmentation function D{z) 



(Ppdy J,, . J d?kdy — . z'^ d?kdy 

where p is the hadron's transverse momentum and 

D{z) oc 1 fl - i - ) , (26) 



z \ z \ — z ^ 

with e = 0.043 for c-quark and e = 0.006 for b-quark. The results of the calculations are 
displayed in Fig. [3] and Fig. HI 




FIG. 3: Nuclear modification factor for open charm production at RHIC and LHC. 

Note, that the calculations of are applicable only if the coherence length of the 

typical qq pair is much larger than the size interaction of the interaction region which, in the 
nucleus rest frame, is about the size of the nucleus for the most central events. It is easy to 
verify (see e. g. 38[) that this condition holds for all lines in Fig. [3] and Fig. H] apart from the 
B-meson production at ?/ = at RHIC. In the former case, the hadron nuclear attenuation 
effect can be responsible for suppression (not show here) of B production. 

In this paper we used a model for A^(r, y) which was tested in the RHIC kinematics. 
Extension to LHC involves an assumption that this model correctly captures the energy 
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FIG. 4: Nuclear modification factor for open beauty production at RHIC and LHC. 

dependence of the inclusive cross sections. It is of great interest to look in yet another 
kinematic region of low x and high m^. The very first paper on geometric scaling 41 1 
hints that the geometric scaling holds in quite a wide kinematic region of a; < 0.01 and 
< 450 GeV^ (and all small x data available). It was explained in |42i] that this property 
stems directly from solution ^BK equation. The geometric scaling is predicted to hold up 
to momenta Qgcom ^ Ql/ ^ [43]. Although there is a fair amount of ambiguity in the value 
of A we can use the model of Sec. [TTlto estimate that at LHC at ?/ = Qgoom = 2.7 GeV, 



and at y = 2 Q^com = 5.0 GeV. On the other hand, the original observation of [4l| as 
well as the parameterization of pj] seems to imply that the geometric scaling may still be 
a reasonable approximation at even higher momenta, provided that x is small enough. In 
order to illustrate dependence of RpA on we calculated it for four quark flavors shown in 
Fig. [H At masses of the order of nit and at rapidity y = 2 at LHC we have x ~ 0.01, so that 
X is indeed small enough for the dipole model to be applicable. The amount of suppression 
of RpA depends on the anomalous dimension 7 as, roughly RpA ~ I/A^^oYm where 7 is a 
function of y and m±, see ([2]) and ([3]). Therefore, by measuring the m±_ and y dependence of 
RpA one can determine the anomalous dimension 7 and map the transition region between 
the geometric scaling and the collinear factorization. The failure of the geometric scaling 
will manifest itself by rapid approach of the RpA to unity (modulo DGLAP corrections) at 
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transverse mass m±_ > Q 



geom • 



Vs=5500 GeV, y=2, k =0.5 GeV 




0.2 - 



10"^ 1 10 10^ 
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FIG. 5: Dependence of the nuclear modification factor on quark mass. Solid line is RpA for quarks 
(no fragmentation). Geometric scaling is expected to break down at m± ~ Qgeom and therefore 
RpA is anticipated to deviate from the solid line towards unity. Dotted lines illustrates a possible 
behavior of RpA- 
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